Healthy LowCD4 AIDS Death
Healthy 0.721 0.202 0.067 0.010
LowCD4 0.000 0.581 0.407 0.012
AIDS 0.000 0.000 0.750 0.250
Death 0.000 0.000 0.000 1.000
Backwards Conversion, Solving for PSA Parameters and Copula-Based PSA Sampling
Transition probability matrix:
Healthy LowCD4 AIDS Death
Healthy 0.721 0.202 0.067 0.010
LowCD4 0.000 0.581 0.407 0.012
AIDS 0.000 0.000 0.750 0.250
Death 0.000 0.000 0.000 1.000
A. Compute the matrix logarithm of the transition probability matrix \(P\).
B. Maximum likelihood-based approach.
C. Convert the supplied transition probabilities to rates one-by-one.
(You can always try 2+ of these methods to see how they compare.)
\[ R = \log P \]
\[A' = V^{-1} A V\]
\[A' = V^{-1} A V\]
\[\log P = V (\log A') V^{-1}\]
V <- eigen(mP)$vectors
iV <- solve(V)
Ap <- iV %*% mP %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
dimnames(R) = dimnames(mP)
round(R,3) Healthy LowCD4 AIDS Death
Healthy -0.327 0.311 0.002 0.013
LowCD4 0.000 -0.543 0.615 -0.072
AIDS 0.000 0.000 -0.288 0.288
Death 0.000 0.000 0.000 0.000
| Healthy | LowCD4 | AIDS | Death | |
|---|---|---|---|---|
| Healthy | -0.327 | 0.311 | 0.002 | 0.013 |
| LowCD4 | 0 | -0.543 | 0.615 | -0.072 |
| AIDS | 0 | 0 | -0.288 | 0.288 |
| Death | 0 | 0 | 0 | 0 |
| Healthy | LowCD4 | AIDS | Death | |
|---|---|---|---|---|
| Healthy | -0.327 | 0.311 | 0.002 | 0.013 |
| LowCD4 | 0 | -0.543 | 0.615 | -0.072 |
| AIDS | 0 | 0 | -0.288 | 0.288 |
| Death | 0 | 0 | 0 | 0 |
expm::logm Higham (2008) method returns same as eigenvalue method.expm::logm() Healthy LowCD4 AIDS Death
Healthy -0.327 0.311 0.002 0.013
LowCD4 0.000 -0.543 0.615 -0.072
AIDS 0.000 0.000 -0.288 0.288
Death 0.000 0.000 0.000 0.000
Steps:
tr0 <- # Initial state occupancy.
c("Healthy" = 1000,"LowCD4" = 0, "AIDS" = 0, "Death" = 0)
tr1 <- # State occupancy after one cycle.
tr0 %*% mP
tr <- # Bind together into a data frame.
rbind.data.frame(tr0, tr1) %>%
mutate(cycle = c(0, 1)) %>%
select(cycle,everything())
tr cycle Healthy LowCD4 AIDS Death
1 0 1000 0 0 0
2 1 721 202 67 10
Healthy state.Healthy while 202 are in LowCD4, etc.tr <-
rbind.data.frame(tr0, tr1) %>%
mutate(cycle = c(0, 1)) %>%
gather(state,count,-cycle) %>%
mutate(count = round(count,0)) %>%
arrange(cycle) %>%
lapply(.,rep,.$count) %>%
cbind.data.frame() %>%
as_tibble() %>%
mutate(state = factor(state,
levels = c("Healthy","LowCD4","AIDS","Death"))) %>%
arrange(cycle,state) %>%
mutate(ptnum = rep(1:1000,2)) %>%
select(-count) %>%
arrange(ptnum,cycle) %>%
mutate(state = as.numeric(state)) %>%
select(ptnum,cycle,state)# A tibble: 2,000 × 3
ptnum cycle state
<int> <dbl> <dbl>
1 1 0 1
2 1 1 1
3 2 0 1
4 2 1 1
5 3 0 1
6 3 1 1
7 4 0 1
8 4 1 1
9 5 0 1
10 5 1 1
# ℹ 1,990 more rows
expm(mP).
Call:
msm(formula = state ~ cycle, subject = ptnum, data = tr, qmatrix = Q.init)
Maximum likelihood estimates
Transition intensities
Baseline
Healthy - Healthy -0.32712 ( -3.680e-01, -2.907e-01)
Healthy - LowCD4 0.32702 ( 1.354e-01, 7.898e-01)
Healthy - Death 0.00010 ( 0.000e+00, Inf)
LowCD4 - LowCD4 -0.64477 ( -1.139e+01, -3.650e-02)
LowCD4 - AIDS 0.63465 ( 9.115e-06, 4.419e+04)
LowCD4 - Death 0.01012 ( 8.238e-299, 1.242e+294)
AIDS - AIDS -0.34837 ( -3.805e+40, -3.189e-42)
AIDS - Death 0.34837 ( 3.189e-42, 3.805e+40)
-2 * log-likelihood: 1572.2
Multistate-Model:
Healthy LowCD4 AIDS Death
[1,] 721 202 67 10
Healthy LowCD4 AIDS Death
[1,] 721 202 67 10
log numerically unstable converting from probability to rate.Multistate-Model:
With a reasonable generator matrix defined, we can now augment the model as we see fit:
Once done with the above, just embed the new matrix.
We have defined our underlying model (either from bottom up or through backwards conversion) but need to define PSA distributions for model parameters.
Model draws on literature-based parameters that are reported as means, standard deviations, interquartile ranges, 95% confidence intervals, etc.
Straightforward to obtain base case values from literature.
But how can we define PSA distributions based on limited information?
| Parameter Type | Distribution |
|---|---|
| Probability | beta |
| Rate | gamma |
| Utility weight | beta |
| Right skew (e.g., cost) | gamma, lognormal |
| Relative risks or hazard ratios | lognormal |
| Odds Ratio | logistic |
ParameterSolver implemented in Windows, Link
See Cook (2010) and Vasco Grilo’s blog post for more.
The next few slides provide you with nearly all the tools you’ll need, however.
For a uniform distribution with minimum \(a\) and maximum \(b\)
\[ a = \frac{p_2x_1 - p_1x_2}{p_2-p_1} \] \[ b = \frac{(1-p_1)x_2-(1-p_2)x_1}{p_2-p_1} \]
\[ \sigma = \frac{x_2 - x_1}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]
\[ \mu = \frac{x_1\Phi^{-1}(p_2)-x_2\Phi^{-1}(p_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]
\[ \sigma = \frac{\ln(x_2) - \ln(x_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]
\[ \mu = \frac{\ln(x_1)\Phi^{-1}(p_2)-\ln(x_2)\Phi^{-1}(p_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]
[1] 299.99
[1] 749.98
[1] 300
[1] 750
\[\min_{\vec{\theta} \in R}\,\, (F(x_1|\vec{\theta})-p_1)^2+(F(x_2|\vec{\theta})-p_2)^2\]
\[\min_{\vec{\theta} \in R}\,\, (F(x_1|\vec{\theta})-p_1)^2+(F(x_2|\vec{\theta})-p_2)^2\]
\[\min_{\vec{\theta} \in R} \sum_{i \in 1,2}\,\, (\text{logit} (F(x_i|\vec{\theta}))-\text{logit} (p_i))^2\]
Tf <- function(x, shape, rate)
pgamma(x,shape,rate,log=TRUE) -
pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)
norm <- function(x1, p1, x2, p2, shape,rate)
(Tf(x1, shape, rate)-(log(p1)-log(1-p1)) )^2 +
(Tf(x2, shape, rate)-(log(p2)-log(1-p2)) )^2
fn <- function(x) norm(x1, p1, x2, p2, x[1], x[2])
gamma_optim <-
optim(c(0.5, 0.1), # initial parameter guesses
fn,
gr = function(x) pracma::grad(fn, x),
lower=c(-Inf, 1e-5),
method = "L-BFGS-B",
control=list(factr=1e-10, maxit=100))$parpbeta rather than pgamma).log=TRUE.In the HIV model discussed earlier, patients accrue two types of costs:
Healthy, LowCD4, AIDS)We can easily draw correlated PSA samples via copulas.
Copulas allow us to sample from the joint distribution of correlated parameters
\[ F(x_1,...,x_d) = C \big ( F_1(x_1), ..., F_d(x_d) \big ) \]